The Hubbard model on the triangular lattice: Spiral order and spin liquid 
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We investigate the half-filled Hubbard model on an isotropic triangular lattice with the variational 
cluster approximation. By decreasing the on-site repulsion U (or equivalently increasing pressure) 
we go from a phase with long range, three-sublattice, spiral magnetic order, to a non-magnetic Mott 
insulating phase - a spin liquid - and then, for U < 6.7t, to a metallic phase. Clusters of sizes 
3, 6 and 15 with open boundary conditions are used in these calculations, and an extrapolation to 
infinite size is argued to lead to a disordered phase at U = 8t, but to a spiral order at U > 12. 



The effect of geometric frustration on quantum mag- 
netism is still a very active field of investigation. The 
quantum Hciscnbcrg model on a two-dimensional square 
(bipartite) lattice exhibits long-range Neel order, but 
that order is suppressed on an isotropic triangular lat- 
tice. In that case, the classical ground state is a spiral 
configuration in which the magnetization on each of the 
three sub-lattices is oriented at 120° of the other two. 
For a while, it was conjectured that quantum fluctua- 
tions around that classical ground state would be strong 
enough to destroy this ordered pattern, but there is now 
a quasi-consensus that this is not the case [1, 2]. The 
latest Monte Carlo studies of the quantum Heisenberg 
model on a triangular lattice point towards a sub-lattice 
magnetization of m w 0.41 in the ground state [2]. 

However, real antiferromagnets are better described by 
the Hubbard model, 

H = -tJ2 C L C J° + Uj2^n ll (1) 

where t is the hopping amplitude between neighboring 
sites, Ciu destroys an electron of spin a at site i and U 
is the on-site Coulomb repulsion. The Heisenberg model 
is recovered in the strong coupling limit (U 3> i), with 
direct-exchange constant J ~ 4t 2 /U. Finite- 1/ effects 
are potentially important on real systems to which the 
Heisenberg model is usually applied. Such effects are 
often incorporated as ring-exchange terms in spin mod- 
els [3], but their origin can be traced back to the Hub- 
bard model itself [4]. For instance, the organic conductor 
k-(BEDT-TTF) 2 Cu 2 (CN) 3 may be described by a Hub- 
bard model on an almost isotropic and half-filled trian- 
gular lattice [5], and this material is conjectured to be 
in a spin liquid (i.e. magnetically disordered and insu- 
lating) phase [6]. So is the triangular antiferromagnet 
EtMe 3 Sb[Pd(dmit) 2 ] 2 [7]. The question that arises in 
this case is whether such a state is compatible with a 
Hubbard model description. In this paper, we will argue 
that it is, i.e., that the Hubbard model on a triangular 
lattice exhibits a spin liquid phase at intermediate values 
of U (e.g. U ~ 8) although it exhibits spiral magnetic 
order at stronger coupling (e.g. at U = 12). 




FIG. 1: (Color online) clusters used in our study. The 6-site 
and 15-site clusters tile the lattice only when paired with iden- 
tical, inverted clusters. Superlattice basis vectors are shown. 

The Hubbard model on an anisotropic triangular lat- 
tice has been studied by various methods. The 120° spiral 
state has been studied in the mean-field approximation 
[8, 9], and a spin stiffness analysis points to a loss of order 
for U < 6 [9], In the isotropic case, slave-bosons methods 
were used to obtain a phase diagram qualitatively similar 
to the Hartree-Fock results [10], with a transition from 
a metallic phase to a magnetic phase, with no interca- 
lated spin liquid phase. On the other hand, the presence 
of a Mott phase was confirmed in Refs. 11-14, however 
without confronting it with a spiral magnetic order. 

In this work we use the (zero temperature) variational 
cluster approximation (VCA) [15]. This method goes be- 
yond mean field and takes into account exactly the effects 
of strong short-range correlations. As U is increased, we 
show that the system goes from a metallic phase to a non 
magnetic, insulating phase (i.e., a spin liquid) at around 
U w 6.7, and then to a magnetic, spiral phase at larger 
values of U. Our treatment involves exact solutions of 
the model on triangular clusters of 3, 6 and 15 sites, as 
well as an extrapolation to infinite size. 

The Variational Cluster Approximation. The VCA 
[16] is a quantum cluster approach to the Hubbard model 
that rests on Potthoff's self-energy functional approach 
(SFA) [15]. It has been applied, for instance, to the prob- 
lem of competing phases in the high-T c cuprates [17] and 
in the layered organic conductors [12]. The SFA involves 
a functional f2t P] of the self-energy, parametrized by the 
one-body terms collectively labeled by t, that is station- 
ary at the physical self-energy of the system: Sfl/6'E = 0. 
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The SFA introduces a reference Hamiltonian H' , with the 
same two-body interaction as the original Hamiltonian 
H, but with a different one-body part, so that H' may 
be solved numerically. The functional f2 t [S] is then 

fitp] =« t 'P]-Trln(Go 1 -S)+Trln(G^ 1 -S) (2) 

where Go and G'o are the non-interacting Green func- 
tions of H and H' , respectively. At the physical self- 
energy, this functional is the grand potential f2. 

In VCA, H' is obtained from H by (i) tiling the lattice 
into a super-lattice of identical clusters, and removing 
all inter-cluster hopping terms and (ii) introducing on 
the clusters Weiss fields that allow for broken symmetry 
phases. Then the Weiss fields (collectively denoted h in 
what follows) are used as variational parameters and the 
functional iltp] reduces to a function flt(h) given by 

nt(t') = n'-/ ^^lndet(l+(G - 1 -G - 1 )G') (3) 

JC k 

where G' is the exact Green function of H' , f2' the exact 
grand potential of H' and the sum is over wave- vectors K 
of the Brillouin zone of the super-lattice. In practice, one 
searches for the stationary points of the above function, 
whose evaluation requires, at each point h, the exact solu- 
tion of the Hamiltonian H' defined on a finite-size cluster. 
At these points, the self-energy £ of H' is considered an 
approximation to the physical self-energy and is used to 
construct the Green function G = (Gg -1 — of the lat- 
tice model. Thus, VCA provides us with an approximate 
Green function of the system, allowing the calculation of 
spectral and thermodynamic properties, both in broken 
symmetry phases and normal phases. 

In investigating broken symmetry phases, VCA is su- 
perior to static mean field approaches in that it does not 
require any factorization of the interaction, and short- 
range correlations (within a cluster) are taken into ac- 
count exactly. The Green function obtained is still de- 
fined on the infinite lattice. The only approximation 
comes from the limited space of self-energies on which 
the variational principle is applied, limited by the cluster 
size and by the number of variational parameters used. 

Clusters for the spiral order. The Weiss field term as- 
sociated with the spiral magnetic order may be expressed 
as H' h — hWl, where 

fBl= J^eA-Si + ^eB-Si+^ec-Si 1 (4) 
[ieA ieB iec ) 

where A, B and G stand for the three sub-lattices of the 
triangular lattice, as shown on Fig. 1 by different shades 
of gray. The unit vectors b.c are oriented at 120° of 
each other, and the spin operator is Sj = cJ cr a( gCi j( g. 

The clusters used in applying VCA to the triangular 
lattice are depicted on Fig. 1. They all have a triangular 
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FIG. 2: (color online) Scaled Potthoff functional £1 as a func- 
tion of Weiss field h for various values of U (3-site cluster). 
The local minima are indicated by arrows. 

shape and treat the three sub-lattices on the same foot- 
ing. Since the L = 6 and L = 15 clusters do not tile the 
lattice by themselves, they are paired with their rotated 
mirror-image to define a true super-lattice, in the Bra- 
vais sense. More explicitly, the Green function G' of the 
super-lattice's unit cell (the union of the cluster and of 
its mirror-image) is given by 

G'- 1 = G' 1 - 1 + G' 2 - 1 +t 12 (5) 

where G[ is the Green function of the cluster itself (site 
and spin indices suppressed), G' 2 that of its rotated mir- 
ror image (a simple transformation of G^) and ti 2 the 
hopping matrix linking the two (dashed links on Fig. 1). 

The variational parameters used in this work are the 
Weiss field h multiplying Wl (see Eq. (4)) and the chemi- 
cal potential // of the cluster. Treating // as a variational 
parameter instead of setting // = /i ensures thermody- 
namic consistency, i.e., that the densities obtained by 
calculating n = TrG and n — —dfl/dfi coincide. Fig. 2 
illustrates the h dependence of Q,t{h,fjf) — ilt(0,[i') for 
the value of // corresponding to the solution, for several 
values of U and on a 3-site cluster. As one can see, a 
local minimum exists as a function of h for U > 5. h and 
fit were divided by J = 4i 2 /U in order to emphasize the 
strong-coupling scaling behavior. 

The Newton-Raphson algorithm is used to locate the 
values of h and // that make the function £lt(h, fi 1 ) 
(Eq. (3)) stationary. The self-energy obtained at that 
point is then used to construct the approximate lattice 
Green function. The spiral order parameter m, i.e., the 
expectation value (SDt) divided by the number of lattice 
sites, is calculated from that Green function as 

m = 2z / j ^ G -»(^,K)Tl ba (6) 

where the frequency integral is taken along the positive 
imaginary axis, and the sum over repeated indices is im- 
plicit. a,b are composite indices including both cluster 
site and spin: a = (i,cr). 9Jl a b is a matrix of real numbers 
expressing Wl as a one-body operator: dJlabC^Cb- 
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FIG. 3: (color online) Left panel: U dependence of the spiral 
order parameter for L — 3, 6 and 15. Right panel: Spiral 
order parameter as a function of scaling parameter Q, for 
various U's. The U — 8 curve is a guide to the eye only. Left 
inset: Neel order parameter as a function of U on a square 
lattice (12-sites). Right inset: The same as a function of Q. 
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FIG. 4: (color online) Scaled Weiss field as a function of Q 
(solid lines) and 1 — 1/1/ (dashed lines) for various values of 
U. The data are obtained for 3- 6- and 15-sites triangular 
clusters. Top left panel: square lattice results at U = 16 for 
the Neel Weiss field, with L = 2, 4, 8, 10, 12 and 16 sites. 



The calculation is performed for several values of the 
lattice chemical potential /i until the density n is close 
enough to half filling (n = 1). In practice, this is easily 
accomplished when U is large enough for a spectral gap 
to open (U ss 6 and above). 

The left panel of Fig. 3 shows the spiral order param- 
eter as a function of interaction strength U/t, for fixed 
cluster size. The order parameter is seen, as expected, 
to saturate at strong coupling. The transition from the 
magnetic to the disordered state seems of first order, but 
the discontinuity depends on cluster size and might dis- 
appear in the thermodynamic limit; we could perform no 
quantitative analysis on this matter. 

The values found for the Weiss field h and the order 
parameter m depend both on U and on cluster size. We 
are naturally interested in the infinite-size extrapolation 
of these, since an ordered solution found on a finite clus- 
ter can disappear in the thermodynamic limit because 
of long wavelength fluctuations of the order parameter. 
Such an extrapolation is very difficult to do with the 
small clusters at our disposal. Moreover, these clusters 
have open, not periodic, boundary conditions. This im- 
plies that the number of sites of the cluster (L) is not the 
only scaling parameter: the size of its boundary could 
also be significant. We define a scaling parameter Q as 
the number of links within the cluster divided by the total 
number of links of the original lattice within a unit-cell 
of the super-lattice of clusters. Q increases with cluster 
size and reaches unity in the thermodynamic limit. It is 
equal to 1/3, 1/2, and 2/3 respectively for the 3-, 6-, and 
15-site triangular clusters. 

We expect the Weiss field to vanish in the thermody- 
namic limit, as it is then no longer necessary to stabilize 
order. If the Weiss field extrapolates to zero at Q < 1, 
this is to be interpreted has a suppression of order due 
to long wavelength fluctuations. Fig. 4 displays the spi- 



ral Weiss field as a function of both Q and 1 — 1/L for 
several values of U, as well as the Neel Weiss field for 
several square lattice clusters at U — 16 and half-filling. 
The square-lattice results show that Q is a better scaling 
parameter than 1 — 1/L since the values of the Neel Weiss 
field neatly fall on a straight line. This line crosses the 
abscissa very close to Q = 1, as it should since long-range 
Neel order is expected in the square-lattice case. In the 
triangular case, however, the 3-site cluster is too small to 
be in the scaling regime, and we must rely only on the 6- 
and 15-site clusters to extrapolate towards Q = 1. Fig. 4 
shows that the Weiss field extrapolates to zero very near 
Q = 1 for all values of U studied except U = 8. Even 
though 1/L scaling looks superficially better for triangu- 
lar clusters, it extrapolates beyond 1/L = (except for 
[7 = 8). Whatever the extrapolation scheme, we con- 
clude that there is no long-range order at [7 = 8. Thus 
long-range spiral order is established somewhere between 
U = 8 and U = 12. 

The right panel of Fig. 3 shows the spiral order pa- 
rameter m as a function of the scaling parameter Q. A 
linear extrapolation to Q — 1 yields m ~ 0.65, which is 
larger than values obtained for the Hciscnberg model by 
Monte Carlo methods [2]. Thus, despite the extrapola- 
tion, VCA still exaggerates the tendency of the system 
to order. Indeed, a similar analysis for the square lattice 
Neel order (inset of Fig. 3) yields an extrapolated mag- 
netization of 0.74, whereas the accepted value is closer 
to 0.61 [18]. Therefore, if VCA predicts a magnetically 
disordered state, it is very likely correct. Incidcntly, it is 
impossible to extrapolate the U = 8 order parameter to 
infinite size, which is a further indication of the absence 
of order at U = 8. 

We now turn our attention to the Mott transition. 
Fig. 5B shows the density of states (DOS), calculated by 



4 




FIG. 5: (color online) Bottom panel (B) : Density of states 
at U = 8 (see text for details). Inset : the Mott gap as a 
function of Q. Top panel (A) : U dependence of the infinite- 
size extrapolated Mott gap. 



integrating the lattice Green function G(u>, K) over wave- 
vectors, for U = 8. The smooth curve is obtained by giv- 
ing the complex frequency an imaginary part 77 = 0.2t, 
equivalent to a Lorentzian broadening of the spectral 
function delta peaks. The red (jagged) curve is the point- 
wise extrapolation towards r\ — > of the DOS calculated 
at r\ = 0.01, 0.005 and 0.002. This extrapolation allows 
for a better numerical estimate of A. This estimate can 
be extrapolated to infinite size (with the help of the scal- 
ing parameter Q, see inset). The extrapolated values can 
be plotted as a function of U to extract a critical value 
U c for the Mott transition (Fig 5A). We find U c » 6.7. 

To conclude, the system is predicted to be a metal for 
U < 6, a magnetically disordered Mott insulator (or spin 
liquid) at intermediate values of U, and a spiral magnet 
at larger values of U (already at U = 12). A coexistence 
of the latter two phases could occur if the transition be- 
tween the two were still of first order in the thermody- 
namic limit, although we cannot conclude on the matter. 
This is in contrast to the square-lattice model, in which 



the Mott transition is pre-empted by Neel order down to 
U = (inset of Fig. 3). 
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